SSA Poverty - Spatial Regression

library(rgdal)
library(foreign)
library(data.table)
library(rhandsontable)
library(spdep)
library(tmap)
# Load districts and attributes
g2 <- readOGR("./data", "svyMaps_2016.06.22_sara")

OGR data source with driver: ESRI Shapefile 
Source: "./data", layer: "svyMaps_2016.06.22_sara"
with 2078 features
It has 11 fields

dt2 <- read.dta("./data/SSApoverty_Dist_forGWR.12.dta")
# Keep STATA labels
dt2.lbl <- data.table(varCode=names(dt2), varLabel=attr(dt2, "var.labels"))
setkey(dt2.lbl, varCode)
dt2 <- data.table(dt2)
# Merge attributes to shapefile
g2.dt <- data.table(g2@data)
g2.dt[, rn22 := as.character(rn)]
g2.dt[, rn := row.names(g2)]
# Recode Ethiopia woredas
dt2[svyCode=="eth2010", svyL2Cd := svyL1Cd * 10000 + svyL2Cd]
setkey(g2.dt, ISO3, svyCode, svyL1Cd, svyL2Cd)
setkey(dt2, ISO3, svyCode, svyL1Cd, svyL2Cd)
# Any unmatched?
dt2[!g2.dt][, .N, by=svyCode]

Empty data.table (0 rows) of 2 cols: svyCode,N

# Remove `rn` from Sara's vars
dt2[, rn := NULL]
# Seems okay, so merge
g2.dt <- dt2[g2.dt]
names(g2.dt)[names(g2.dt) %like% "rn"]

[1] "rn"   "rn22"

# Neighbor list
coords <- coordinates(g2)
nb2 <- poly2nb(g2, row.names=paste(g2$ISO3, g2$rn, sep="."))
nb2

Neighbour list object:
Number of regions: 2078 
Number of nonzero links: 10902 
Percentage nonzero weights: 0.2524731 
Average number of links: 5.246391 
3 regions with no links:
MOZ.424 MWI.757 MWI.780

Plot contiguities for e.g. Ghana

# Plot contiguities (e.g. GHA)
gha <- g2[g2$ISO3=="GHA",]
gha.dt <- g2.dt[ISO3=="GHA"]
gha <- SpatialPolygonsDataFrame(gha, data.frame(gha.dt), match.ID="rn")
coords.gha <- coordinates(gha)
nb2.gha <- poly2nb(gha)
plot(gha, border="grey",
  main=list("Contiguity - Ghana (170 districts)"))
plot(nb2.gha, coords.gha, col="blue", add=T)

legend("bottomleft", 
  legend=capture.output(summary(nb2.gha))[1:7], cex=.6)

# Spatial weights for SSA (and Ghana only)
#w2 <- nb2listw(nb2)
w.gha <- nb2listw(nb2.gha)

# Explanatory vars
var <- c("spei_lt", "L1_speihishock", "L1_speiloshock")

# Plot to identify districts with missing values (Ghana only)
tmap_mode("view")
qtm(gha, fill=var[1])

# Impute missing values with regional median
tmp <- gha.dt[, lapply(.SD, median, na.rm=T), by=svyL1Cd, .SDcols=var]
setnames(tmp, 2:4, paste(names(tmp), "imp", sep="_"))
setkey(tmp, svyL1Cd)
setkey(gha.dt, svyL1Cd)
gha.dt <- tmp[gha.dt]

gha.dt[is.na(spei_lt), spei_lt := spei_lt_imp]
gha.dt[is.na(L1_speihishock), L1_speihishock := L1_speihishock_imp]
gha.dt[is.na(L1_speiloshock), L1_speiloshock := L1_speiloshock_imp]

# Spatial weights for SSA (and Ghana only)
#w2 <- nb2listw(nb2)
w.gha <- nb2listw(nb2.gha)
# Explanatory vars
var <- c("spei_lt", "L1_speihishock", "L1_speiloshock")
# Plot to identify districts with missing values (Ghana only)
tmap_mode("view")

tmap mode set to interactive viewing

qtm(gha, fill=var[1], colorNA="grey", textNA="Missing")

# Impute missing values with regional median
tmp <- gha.dt[, lapply(.SD, median, na.rm=T), by=svyL1Cd, .SDcols=var]
setnames(tmp, paste(names(tmp), "imp", sep="_"))
setkey(tmp, svyL1Cd_imp)
setkey(gha.dt, svyL1Cd)
gha.dt <- tmp[gha.dt]
gha.dt[is.na(spei_lt), spei_lt := spei_lt_imp]
gha.dt[is.na(L1_speihishock), L1_speihishock := L1_speihishock_imp]
gha.dt[is.na(L1_speiloshock), L1_speiloshock := L1_speiloshock_imp]
gha <- SpatialPolygonsDataFrame(gha, data.frame(gha.dt), match.ID="rn")
# Moran plots of explanatory vars (e.g. Ghana)
for(i in var) {
  moran.plot(gha@data[[i]], w.gha, zero.policy=T,
    xlab=dt2.lbl[i, varLabel], 
    ylab=paste("Spatially Lagged", dt2.lbl[i, varLabel]),
    labels=gha$svyL2Nm
  )
}

# Compare models
m1 <- lm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt)
m1s <- sacsarlm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt, w.gha, zero.policy=T)

# Results
summary(m1)
summary(m1s)

# Label 20 districts at random
rnd <- sample(1:170, 20)

# Compare models
m1 <- lm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt)
m1s <- sacsarlm(pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock, 
  data=gha.dt, w.gha, zero.policy=T)
# Results
summary(m1)

Call:
lm(formula = pcexp_ppp_m ~ spei_lt + L1_speihishock + L1_speiloshock, 
    data = gha.dt)

Residuals:
    Min      1Q  Median      3Q     Max 
-118.41  -38.20   -9.39   24.87  288.70 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     165.165      5.291  31.217  < 2e-16 ***
spei_lt         246.757     53.750   4.591 8.92e-06 ***
L1_speihishock  -80.199     37.767  -2.124   0.0353 *  
L1_speiloshock   -8.735     37.650  -0.232   0.8168    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 64.21 on 159 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1306,    Adjusted R-squared:  0.1142 
F-statistic: 7.964 on 3 and 159 DF,  p-value: 5.573e-05

summary(m1s)

Call:
sacsarlm(formula = pcexp_ppp_m ~ spei_lt + L1_speihishock + L1_speiloshock, 
    data = gha.dt, listw = w.gha, zero.policy = T)

Residuals:
      Min        1Q    Median        3Q       Max 
-116.6812  -37.7663   -9.7939   24.8756  288.1350 

Type: sac 
Coefficients: (asymptotic standard errors) 
               Estimate Std. Error z value  Pr(>|z|)
(Intercept)    138.1537    55.0222  2.5109   0.01204
spei_lt        235.3452    57.4579  4.0960 4.204e-05
L1_speihishock -81.2082    36.3581 -2.2336   0.02551
L1_speiloshock  -7.2666    37.0376 -0.1962   0.84446

Rho: 0.16848
Asymptotic standard error: 0.34074
    z-value: 0.49446, p-value: 0.62098
Lambda: -0.16499
Asymptotic standard error: 0.39596
    z-value: -0.41669, p-value: 0.67691

LR test value: 0.14552, p-value: 0.92982

Log likelihood: -907.6312 for sac model
ML residual variance (sigma squared): 3972.9, (sigma: 63.031)
Number of observations: 163 
Number of parameters estimated: 7 
AIC: 1829.3, (AIC for lm: 1825.4)

# Label 20 districts at random
rnd <- sample(1:170, 20)
# simple LM
plot(gha.dt[!is.na(pcexp_ppp_m), pcexp_ppp_m], m1$fitted.values,
  main="Simple LM", 
  xlab="pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock",
  cex=.5, pch=16)
text(gha.dt[!is.na(pcexp_ppp_m)][rnd, pcexp_ppp_m], m1$fitted.values[rnd], 
  labels=gha.dt[!is.na(pcexp_ppp_m)][rnd, svyL2Nm],
  cex=.6, pos=4)

# SARAR
plot(m1s$y, m1s$fitted.values,
  main="SAC/SARAR LM", 
  xlab="pcexp_ppp_m~spei_lt+L1_speihishock+L1_speiloshock",
  cex=.5, pch=16)
text(m1s$y[rnd], m1s$fitted.values[rnd], 
  labels=gha.dt[!is.na(pcexp_ppp_m)][rnd, svyL2Nm],
  cex=.6, pos=4)